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ABSTRACT 

We study the evolution of bound pairs of star clusters by means of direct N-body 
simulations. Our simulations include mass loss by stellar evolution. The initial condi- 
tions are selected to mimic the observed binary star cluster NGC 2136 and NGC 2137 
in the Large Magellanic Cloud. Based on the rather old ages (~ 100 Myr), masses, 
sizes of the two clusters and their projected separation, we conclude that the clus- 
ter pair must have been born with an initial separation of 15-20 pc. Clusters with a 
smaller initial separation tend to merge in ^ 60 Myr due to loss of angular momen- 
tum from escaping stars. Clusters with a larger initial separation tend to become even 
more widely separated due to mass loss from the evolving stellar populations. The 
early orbital evolution of a binary cluster is governed by mass loss from the evolving 
stellar population and by loss of angular momentum from escaping stars. Mass loss 
by stellar winds and supernovae explosions in the first ^ 30 Myr causes the binary to 
expand and the orbit to become eccentric. The initially less massive cluster expands 
more quickly than the binary separation increases, and is therefore bound to initiate 
mass transfer to the more massive cluster. This process is quite contrary to stellar 
binaries in which the more massive star tends to initiate mass transfer. Since mass 
transfer proceeds on a thermal timescale from the less massive to the more massive 
cluster, this semi-detached phase is quite stable, even in an eccentric orbit until the 
orbital separation reaches the gyration radius of the two clusters, at which point both 
clusters merge to one. 

Key words: gravitation - stellar dynamics - methods: iV-bodysimulations - galaxies: 
star clusters - globular clusters: individual: NGC2136 - globular clusters: individual: 
NGC2137 



1 INTRODUCTION 

The large Magellanic cloud (LMC) is a rich environment 
with m any relatively young star clusters iMackev fc Gilmord 
Several of these clusters appear closer together on 
the sky than expected fro m statistical arguments, i.e, they 
look like binary clusters jBhatia et alj Il99lt iDieball et alJ 
12002). The LMC contains a total of 69 cluster pairs with a 
separation R ^ 1'.3 iBhatia fc Hatzidimitrioulll988h which 
corresponds to 19 pc by adopting distance modulu s 



of 18.5 mag to the LMC iGroene'wreg en_fc_Saljji£ 



many of which have similar colors jKontiza^^T 
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The small Magellanic clou d also contains a relatively larg e 
number of paired clusters iHatzidimitriou fc Bhatialll990H . 
Binary star clusters ar e expected to merge on a time scale 
of a few times 10 Myr iBha tia fc Hatzidimitrioulll98d) . and 
NGC 2214 may be an LMC cluster in the process of merging 
('Bha tia fc MacGillivrav 198a) (but see .Banks et al.. (.19951) 
for counter arguments). 

Well known examples of such binary clust ers are 

SL 538/NG C 2006 (SL 537) JPieball fc Grebell 1199^1 . 
NGC 1850 iFischer et alJ 1199311 and NGC 213 6 (SL762, at 
a = 5^ 53" 17"; S = -69° 32")/NGC2137 iHilker et alJ 
^9^. Confirmation of the binarity of such apparent pairs 
of star clusters should come from detailed measurements of 
their orbital parameters, rather than statistical arguments. 
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However, with an orbital period of tens of million years 
such measurements are not trivial. In th e Antennae galaxies 
l|Whitmore et al.l 12005*: 'Fall ct al."2005'^ and in the youn g 
starburst galaxy M51 (Larson 2000: Bastian et alJ boOSl) . 
star clusters appear to be formed in groups often consist- 
ing more than two. Apparently, star clusters in general are 
born in conglomerates in which several clumps tend to form 
within a relatively small time interval presumably initiated 
by galaxy-galaxy interactions. The individual star clusters 
in such conglomerates may subsequently merge in due time 
jFellhauer fc Kroupalliooa iFellhauer et alj|2002ll . Cluster 
binaries may form in a comparable environments, in which 
case their ages have to be at least comparable. Alternatively, 
a more massive cluster could possibly captures a lower mass 
cluster in a tidal event, resulting in two cluster with differ- 
ent ages. However, this process is likely to result in a single 
cluster rather than a binary (jje Olivcira ct al. 1999). 

In this study we mainly concentrate on the two star 
clusters NGC 2136 and NGC 2137, which attract attention 
because the p ro.iected distance between them is only about 
20 pc (1'.34) iStein et alJll994^ . In addition, both clusters 
have similar ages of abo ut 100 ± 20 Myr j birsch ct al. 2000) 
and identical metalicity jHilker et al.ll995F^^ With a p rimary 
mass of 26300-28200 jMackev fc Gilmor'3l2003al) and a 
mass ratio q ~ 0.17, their orbital period is about 46 Myr. 
If bound, both clusters have orbited each-other about twice 
within their lifetime. Regretfully it is currently not known 
whether or not the clusters really form a bound pair, but at 
the moment this seems to be the most logical conclusion. 

In this study we explore the initial conditions for 
which a pair of star clusters could survive until an age 
of about 100 Myr. We focus in particular on the observed 
pair NGC 2136 and NGC 2137 jHilker et al.lll995l) . as they 
form an observational motivation for this exercise. Assum- 
ing that the two clusters form a bound pair we are able to 
limit the initial cluster parameters. We identify two main 
regimes in binary cluster evolution. Small initial orbital sep- 
aration ( ^ 12 pc) tends to lead to a merger on a time scale 
of less than about 60 Myr, whereas large initial separation 
( ^ 17 pc) causes the two star clusters to recede. Clusters 
with initially a larger separation become very vulnerable 
to disruption by passing other star clusters, giant molecu- 
lar clounds or the background potential of the host galaxy. 
Between about 12 pc and 17 pc is a range of orbital separa- 
tions between which the period changes little and in which 
the cluster binary is able to survive for more than 100 Myr. 
In these cases the change in orbital angular momentum by 
stellar mass loss and escapers is compensated by the redis- 
tribution of internal angular momentum in the rotation of 
the clusters. 



2 VALIDATION AND ANALYTIC 
CONSIDERATIONS 

We perform direct A^-body simulations of binary star clus- 
ters using the kira integrator of the starlab simulation 
environment (Portegies Zwart ct al. 2001). This A''-body 
code computes inter-particle forces by direct summation 
and the integration of the equations of motion is carried 
out with a fourth-order H ermite predictor-correct or scheme 
jMakino fc Aarsethlll992h with block time steps dMcMillanl 



Il986h . The greatest spee d is obtained with the GRAPE-6 
speci a l purpose compu ter jMakino et alll997ll20o3 : lMakind 
120041: iMakino fc Taiii Il99d) and we use the MoDeStA^ 
platform in Amsterdam. Some of our calculations incor- 
porate the ev olution of single stars and binar ies via the 
SeBa package jPorteeies Zwart fc Verbundll99^ . However, 
since our simulations all started with single stars, bina- 
ries play a minor role as only a few are formed during 
the course of the simula tions by 3-body dynamical capture 
jAarseth fc Heggidll97dl . 

Following the stellar evolution package SeBa, stars 
more massive than about 25 Mq turn into black holes in 
a supernova explosion, and stars more than about 8 Mq 
turn into neutron stars. The latter receive a kick veloc- 
ity upon formation from the Paczin sky-Hartman d istribu- 
tion with a dispersion of 300 km s^ "^ llHartmanlll997r) . black 
holes of mass mbh receive kicks from the same distribution 
but with a dispersions of 300 x 1 .4M(7)/mbhkm s~^ (see 
IPortegies Zwart fc YungelsonI il998l) for details). During the 
simulation we keep track o f the individual clusters by using 
a clump-finding algorithm iEisenstein fc Hutlll998ft . Masses 
referring to the individual clusters when discussing the re- 
sults in Tab.|5|are therefore an over estimate as the potential 
of a possible companion cluster is ignored in this procedure. 

2.1 Binary-cluster evolution without stellar 
evolution 

To test the merger time for bin ary clusters we started by 
reproducing t he calculations of llSugimoto fc Making. 19891 : 
iMakino et alJ 119911. They performed for t hat time large 
scale N-body simulations using a tree code dBarnes fc Huti 
Il986l) with a softening of e = 0.025rvir. Here rvir is the ini- 
tial virial radius of the primary cluster. All stars in their 
models had the same mass and they ignored stellar evo- 
lution . The initial density profile for all simulations was a 
iKind (Il966l) model with Wo = 7, the tidal radius is then 
rtide — 7rvir. When on a circular orbit with separation 
a = 10 Tvir both clusters are initially in Roche-lobe contact. 
We summarize the initial conditions and results in Tab.0 
These simulations are carried o ut in dimension less A''-body 
units jHeggie fc Mathieulll986h in which M = G = rvir = 1. 
Here M is the total cluster mass and G is Newton's constant. 

When adopting the sa me initial conditio ns as 
ISugimoto fc Makind (Il989f) and IMakino et alJ (Il99lll we re- 
produced their results. However, when adopting zero soften- 
ing £ = our clusters tended to merge somewhat earlier and 
the merger product was somewhat less massive (see TablH 
for details). This is not unexpected as softening in the sim- 
ulations has the effect of reducing close encounters causing 
the relaxation time for the cluster pair to increase. As a re- 
sult, softening tends to reduce mass loss, by reducing the 
number of close encounters, and increases the merger time. 

2.2 Analytical considerations 

The relation between mass loss and angular momentum loss 
from the binary can be computed from the conservative as- 
sumption that both clusters are in synchroneous rotation. 

^ see|http : / /modesta ■ science ■ uva ■ nl | 
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Table 1. Overview of the initial conditio ns explored by 
LSugimoto & Makino 1 1989); Makino et alJ il991^ and reproduced 
here. In the first column is the name of the simulation followed 
by the number of particles in two clusters. The subsequent two 
columns give the initial separation with which the clusters were 
born {Ri) and their initial orbital eccentricity (e;). Note that 
cluster B was born near semi-latus rectum. The last four columns 
give the results of the sim ulations. These give the t ime of merger 
in standard A'^-body units <Heggie fc Mathiejl98^ and the 
total number of stars in the merger cluster as fraction of the total 
number of stars in the initial cluster /jy = [Nj + nf)/{Ni + n;). 
The last two columns give these numbers for the simulations with- 
out softening. 



Model 


TV 


n 


Ri 




*/ 


fN 


*/ 


In 












e = 


0.025rvir 


e 


= 


A 


2048 


2048 


10 


0.0 


370 


0.95 


270 


0.93 


B 


2048 


2048 


6 


0.5 


190 


0.96 


150 


0.95 


C 


4096 


2048 


6 


0.0 


105 


0.97 


100 


0.96 


D 


4096 


2048 


6 


0.0 


145 


0.97 


120 


0.96 



J — fia UJ : 



a ujM. 



+ 



(1) 



Here /x = mpms/M , with M = nip + rris and q = rris/mp. 
The full differential can then be written as 

D\nJ = dln^i + 2dlna + dlnuj. (2) 

With Kepler's iKepieJ dieOSTl third law (a^ oc M/w^) 

3dlna = dlnA/ - 2dlntj, (3) 

we cen reduce Eq.|2|by eliminating the orbital velocity 

din J = din nip + din nis — ^dln M + ^dln a. (4) 

With the assumption that mass is lost from the binary with 
angular momentum 

dJ/dM = jJ/M (5) 
Eq.|3 reduces to 

dlna = (27 + l)dlnM - 2dlnmp - 2dlnms. (6) 

Which, after integrating from the initial conditions to the 
final conditions results in 

(7) 



ai \Mi) 



nip ilTls i 



If mass is predominantly lost by stellar evolution and if 
both clusters have the same initial mass function, we may 
adopt that the mass ratio q remains constant throughout 
the mass loss process. In this case Eq.|21 reduces to 



d In J = d In M + 2d In a + d In tj, 



(8) 



which after making the same assumption about the amount 
of angular momentum lost per unit mass reduces to 



Integrating between the initial and final parameters then 
gives 



Of 

a, 



(Ml)' 



(10) 



Since in our simulations the clusters are initially not ro- 
tating we then overestimate the rate of angular momentum 
loss. 

The angular momentum of two point masses in a circu- 
lar orbit is 



If mass is lost adiabatically and isotropically (i.e, 7=1) 
Eg . 1101 reduces to the classical result in which aM — con- 
stant. 



2.2.1 Results for the analytic prescription 

With the above prescription we can calculate the orbital 
evolution of a binary star clusters in the absence of dynam- 
ical friction. In the presence of stellar evolution most mass 
is likely to be lost via the evolving stellar population. 

To calculate the evolution of the cluster mass we adopte 
an initial mass function, for which we use a lSalpeteJ lll955l) 
between a minimum mass of m_ and a maximum of m+. 
We now assume that all stars with a mass m* > mto are 
lost from the clusters. The lifetime of a star of mass m* is 
cal culated using the fit to detailed stellar evolution tracks 
by jEggleton et alJll989ft 



t*(m*) 



3600 + MOml-^ + 1.4mt^ 
0.033mi-5 + 0.35m^-^ 



(11) 



dlna = (27 - 3)dlnM. 



(9) 



According to this methodology stellar remnants are lost 
from the cluster. As long as the cluster is younger than about 
40 Myr this assumption does not result in an overestimate of 
the mass loss since black holes and neutron stars are likely 
to receive a velocity kick upon formation. As a result the 
majority of compact objects are ejected from the cluster in 
an early stage. At later time, when white dwarfs form, the 
use of Eg . llll tends to overestiemate the mass loss, as white 
dwarfs tend to stay in the cluster. In fig.^we present the 
evolution of the cluster mass for several choises of m_. 

In fig.|21we present the evolution of the orbital separa- 
tion under several assumptions for the amount of angular 
momentum carried away per unit mass (see Eq.|7|l. These 
calculations are performed under the assumption that the 
time scale for mass loss is long compared to the orbital pe- 
riod. We will later see that this is not a valid assumption in 
the first few tens of million years. 



2.2.2 Results for two-body integration 

To test the above derivation we implement the mass loss 
prescription adopted in 5 12. 2. H and integrate the equations of 
motion of a two-body system. During integrating the equa- 
tions of motion we allow the two point masses to lose mass. 
Note than since the prescription for stellar mass loss depends 
on physical scales we provide a scaling to stellar units to the 
2-body integration. 

To enable a direct comparison with the simulations re- 
ported in §|3|we perform our calculations with similar ini- 
tial conditions. In figsj^and 2] we present the evolution of 
the orbital parameters {a and e, respectively) for two point 
masses. The initial orbit was circular, and we adopted a pri- 
mary mass of 27500 M0 and a mass ratio of g = 1.0 and 
g = 0.16. 

The mass loss prescripti o n desc ribed in i) l2.2.1l is then 
combined with the ISalpeteil ^1955^ initial mass function. 
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Figure 1. The mass of a star cluster as a function of time. From 
top to bottom the varous curves are calculated with a Salpeter 
initial mass function with a maximum mass of m+ = 100 Mq 
and with a minimum mass of m_ = 0.04 Mq, 0.2, 1.0, 5 and 
m_ = 25 Mq. We adopted Eg . 1111 to calculate the lifetime of a 
star with mass (see SI2.2.1l for details). 
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Figure 3. Semi-major axis vs. time for various initial conditions 
of two mass losing particles in an initialy circular orbit. The top 
(4) curves are calculated with m_ = 1.0 Mq and the bottom 
curves with m_ = 0.5 Mq. The thick curves are for = 20 pc, 
the thin curves for ai = 10 pc. The solid curves are for q = I, and 
the dotted curves for 5 = 0.1. 
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Figure 2. The variation in the orbital separation as a function of 
time for various choises of 7 = 0, 1, 1.5 and 7 = 2 (see Eq.ljJ. Each 
calculations is performed with a Salpeter initial mass function 
between m_ = 1 Mq and m+ = 100 Mq. 



from which we then can calculate the rate of mass loss. Cal- 
culations are performed with a maximum mass of m= — 
100 Mq and with a minimum of m_ = 1.0 Mq and m_ = 
0.5 Mq. The evolution of the cluster mass for these param- 
eters are presented in Fig.0 

The orbital expansion (Fig.j^J is comparable to that of 
the analytic prescription using 7 = 1 (see Fig.|2l . The semi- 
major axis for the choise of m_ = 1 Mq increases more 
quickly than for m_ = 0.5 Mq. In this two-body approach, 
the friction between the clusters and possible transfer of 
stars from one to the other are ignored. As a result, the 
initial mass ratio and the initial orbital separation have little 
influence on the lifetime of the binary cluster. For orbital 
separations smaller than about 10 pc, it will turn out that 
dynamical friction has an important effect on the orbital 
evolution. We will quantify this in §|21where we perform the 
full A'^-body simulation of the binary star cluster. 

In Fig.^Jwe show the evolution of the orbital eccentric- 




Figure 4. Evolution of the orbital eccentricity for the two or- 
biting clusters with various initial conditions of two mass losing 
particles in an initialy circular orbit (see S I2.2.2I . The top (4) 
curves are calculated with ai = 20 pc and the bottom curves with 
ai = 10 pc. The thick curves are for m_ = 1.0 Mq and the thin 
curves for m_ = 0.5 Mq. The solid curves are for <? = 1, and the 
dotted curves for q = 0.1. 



ity of two orbiting objects. The initially circular orbit picks 
up a considerable eccentricity due to the copious mass loss 
in the first 10 Myr. 

At later time, for several simulations around 20 Myr, 
the orbital eccentricity reduces again. This is a conseqence of 
mass loss in an eccentric orbit. If mass is lost near apocenter 
the orbit tends to become less eccentric whereas mass loss at 
pericenter tends to increase the eccentricity. An orbit with 
a smaller initial mass ratio tends to become somewhat more 
eccentric than for a large mass ratio. 

Our assumption that mass is lost adiabatically, clearly, 
is not supported by the 2-body simulations (see fig.^J- Ac- 
cording to that hypothesis the orbit should remain circular. 
But mass in the first few 10 Myr is lost at such a high rate 
that it should in part be treated as an impulsive event. 
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Impulsive mass loss is a two body system tends to re- 
duce the absolute value of the potential energy and the or- 
bital kinetic energy by decreasing the mass of the system. 
Such sudden mass loss imposes no torgue on the system, 
so the angular momentum per unit mass remains the same, 
but the total angular momentum is reduces as a result of 
the decrease in the mass. Following ISills 1983) we can then 
calculate the orbital eccentricity if the total mass of the bi- 
nary cluster drops from Mi to Mf. Assuming that the initial 
orbit is ciruclar = the final eccentricity becomes llHillsl 
Il983l) . 

_ 1-M,/M. 

M,/M, ■ ^ ' 



3 THE iV-BODY SIMULATIONS 

We now focus on the A'^-body simulations of binary star 
clusters. With the analytic analysis in S I2.2.1I and by inte- 
grating the 2-body systems in S lTT^ we have acquired some 
qualitative understanding of the range of initial conditions 
to search in order to reproduce the observed cluster pair 
NGC2136/NGC2137. 

3.1 Setting up the realistic simulations 

The main problem in generating proper initial conditions for 
NGC 2136/NGC 2137 is that we only approximately know 
the conditions at an age of about 100 Myr. In our simulations 
we will attempt to start with reasonable initial conditions 
and see what effect small changes to those have on the final 
configuration, and whether or not the final system has any 
resemblance with the observed cluster pair. We explain here 
how to obtain such initial conditions. 

Assuming an initial mass function (Salpeter), age (80- 
120 Myr) and metallicity the mass of the primary cluster 
(NGC 2136) can be estimated from its observed luminos- 
ity (B = 10.99, V = 10.7) and turns out to be about 
M = 26300-28200 Mq jMackev fc GilmordEooiah . The sec- 
ondary cluster (NGC 2137) is then with B = 12.66, V = 
12.88 about m = 4500 Mq, i.e. the mass ratio q = m/M ~ 
0.167. At a distance of about 50kpc and with a separation 
between the two clusters of 1'.34, the current projected sep- 
aration between the two cluster is ii ~ 20 pc. Assuming that 
their orbit is circular we can compute the tidal radii of the 
two clusters by estimating the size of their respective Roche 
radii with the approximated equation of EaElcton (198.|). 
With the adopted parameters this result in a tidal radius for 
the primary cluster of Rt ~ 15.0 pc, and for the secondary 
Vt — 5.0 pc. T he observed core radius of the primary cluster 
Rc ~ 2.0 pc jMackev fc GilmorelEo03a^ . It turns out that 
if we assume that the tidal radius for the primary cluster 
(NGC 2136) equals to its Roche radius, it has an unusu- 
ally shallow density pr ofile which cannot be described by a 
iKind lll966ri model or a lPlummeil lligilT) sphere. During the 
cluster lifetime its structure is likely to have changed quite 
substantially and we are unable to determine the initial den- 
sity profile. For consistency with S I2.1l we adopt non-rotating 
King Wo ~ 7 density profiles as the initial conditions for 
each of the two clusters in a binary. 

We set up the simulations by determining the orbit of 



the two clusters. For clarity we adopt here model A_20 as an 
example, initial conditions for the other simulations are pre- 
sented in Tab. 121 We start by generating parameters for two 
point masses {Mi = 27500 M© and Mi = 4400 Mq , for the 
primary and secondary cluster, respectively) in a circular or- 
bit with a orbital separation of Ri = 20 pc. Later we replace 
the point masses with the two clusters. Each of the clus- 
ters is generated by selecting a number of stars A''^ = 9000 
for the primary and rii = 1500 for the secondary cluster, 
both of which are sprinkled in a King (1966) model with 
Wq = 7 and a virial radius of r^ii — 2.14pc for the primary 
and rvir — 0.71 pc for the secondary cluster. With these 
parameters both clusters are initially precisely filling their 
respective Roche lobes. We subsequently assign a random 
mass to each of the stars from a Salpeter initial mass func- 
tion between m_ — IMq and 100 Mq, leading to a total 
mass of about M ~ 27500 Mq for the primary cluster and 
m ~ 4400 Mq for the secondary. We call this models A_20. 
Additional simulations are performed with larger total ini- 
tial mass to compensate the mass loss from the evolving stel- 
lar populations. These other models (called B, C and D) are 
computed with a primary cluster mass of Mi ~ 40800 and 
rrii ~ 6800. We repeated these simulations with Ri vary- 
ing from 10 to 20 pc, and we vary the minimum mass of 
the initial mass function from m_ — 1.0 M© (model B), 
m_ = 0.63 Mq (model C) to m_ = 0.50 Mq (model D). To 
study the effect of the initial mass ratio we performed sev- 
eral additional simulations with q — 0.3, q — 0.6 and q — 1. 
The initial conditions are summarized in Tab.|5| 

3.2 The evolution of simulation models A 

We will now discuss the results for model A, and continue 
discussing models B, C and D in S I3.3I The mass lost by 
the primary and secondary clusters in model A_20 and D_20 
are presented in fig.El For model C_20 we only present the 
result for the primary cluster. Over-plotted with thin curves 
are the results of a semi-analytic calculation presented in 
fig-Q The two thin curves in fig.|Hl follow the evolution of 
simulation models A_20 (and C_20) quite satisfactorily. 

Mass loss from both clusters in the first few million years 
is negligible as only the most massive stars lose mass on the 
main sequence. After that the cluster loses mass rapidly, and 
its effect becomes quit important. By about 20 Myr the clus- 
ters have lost about 20% of their mass. The cluster mass loss 
closely follows the semi-analytic description we presented in 
Apparently little mass is lost in the form of escaping 
stars, consistent with the earlier simulations presented in 

As we discussed in S I2.2I the copious mass loss in the 
first few tens of million years drives an expansion of the 
orbit and, since the time scale for mass loss and the orbital 
time scale are comparable the orbital eccentricity increases. 

In simulation A_20, for example the combined effect of 
the adiabatic/impulsive mass loss induces an eccentricity of 
about e ~ 0.25 in the first ~ 20 Myr, and the fact that the 
cluster approaches apocenter causes the distance between 
the two clusters to increase from 20 pc initially to ~ 31 pc 
at t ~ 20 Myr. The induced eccentricity on the orbit is con- 
sistent with Eg . 1121 1' see also fig.QJ, and also the evolution of 
the orbital separation is consistent with the 2-body simula- 
tions in § 12.2.21 
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Table 2. Initial conditions for the simulations. The first column gives the model name, followed by the initial conditions; initial orbital 
separation (Ri), minimum mass of the Salpeter mass function number of stars in the primary (Ni) and secondary clusters {rii) 

and the virial radii of the primary (Hvir) and secondary clusters (rvir). The next set of columns give the final parameters, time at which 
we stopped the simulation (i/), the number of stars belonging to the primary cluster A^^ and secondary cluster rif and the distance 
between the clusters Rf. Since none of the A_10 models survive the last two columns give the moment of merger and the number of 
stars in the merger. All simulations ware performed with King models Wq = 7. The simulation naming is chosen as follows: model A, 
B, C and D identify the minimum mass for the initial mass function, followed by the initial distance in parsec. All models are computed 
with a mass ratio of g = 0.167, except if noted otherwise in the model name, in which case the number after the Q identifies the adopted 
initial mass ratio. 
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Figure 5. Evolution of the mass for the primary and secondary 
clusters in models A_20 (thick solid curves) and for models D_20 
(thick dashes) and for the primary cluster in simulation model 
C_20 (thick dash-3-dots). The error bars to the right indicate the 
best e estimate for the m asses of clusters NGC2136 and NGC2137 
from lDirsch et ahl i200(J) . Over-plotted with a thin solid and thin 
dash-3-dotted line are the result of our semi-analytical calculation 
using two point masses for model A_20 and C_20. The mass loss 
rate (from S I2.2.1I matches the full N-body simulations quite well 
and the resulting orbital variation is presented in fig. llUI 



In fig.[7|we present a representation of the stellar posi- 
tions for simulation model A_10 (top panels) and A_20 (bot- 
tom) at various moments in time. We stop both simulations 
when they reach an age of at least ~ 80 Myr. The three 
cluster pairs of model A_10 have merged before that time, 




t [Myr] 

Figure 6. Evolution of the separation for simulation models 
A with various initial orbital separations (solid curves). The two 
curves for model A_20 are for the simulation model without stellar 
mass loss (dotted curve) and with stellar mass loss. All other 
curves include stellar mass-loss, but start at a different initial 
separation, varying between 10 pc and 20 pc. The bullets in the 
top and bottom solid curves indicate the moments in time at 
which a snapshot is presented in Fig.lTI 



whereas the distance between the two clusters of model A_20 
has then increased beyond 50 pc. 

The simulations in which the clusters are initially rather 
widely separated can satisfactorily be described by the semi- 
analytic analysis in S I2.2I In close proximity the transport of 
angular momentum from the orbit to the internal rotation 
becomes gradually more important (see i) l2.H : tighter clus- 
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Figure 7. Representation in the X-Y-plane for models A_10 (top panels) and A_20 (bottom panels). The simulations models A_10 
are given at birth (left image), at an age of about 15Myr, 30Myr, 45 Myr and at 60Myr (right), whereas for simulation model A_20 
the images represent the cluster birth (left image), at an age of 21 Myr, 42 Myr, 63 Myr and at 84 Myr (right). The moments at which 
snapshots were taken are identified with bullets in Fig.|H] {Note concerning this LANL version: the published version will contain high 
guality images.) 



ter binaries are able to transport angular momentum more 
effectively from the orbit to the rotation of the individual 
clusters, but also to steUar escapers. These dynamical inter- 
ferences dramatically reduces the predictability of a simple 
semi-analytic model, because it is not apriory clear what 
amount of specific angular momentum in stored in the clus- 
ter. For close pairs, stellar mass loss appears to be relatively 
less important whereas the redistribution and loss of angu- 
lar momentum dominates the evolution of the orbital pa- 
rameters. This is illustrated in Fig.jS] where we present the 
evolution of the orbital separation for simulation models A 
with an initial separation of 20 pc, 17.5 pc, 15 pc, 12.5 pc and 
10 pc (solid curves from top to bottom). The orbits of the 
clusters at 10 pc, 12.5 pc and 15 pc expand less dramatically 
than expected on the amount of mass lost from the evolving 
stellar populations. 

The separation between the two clusters in models A_15 
increases in the first ~ 40 Myr by mass loss and the orbital 
eccentricity increases but the loss of angular momentum pre- 
vent that the clusters continue to recede, as is the case in 
model A_20. In the latter model the two clusters approach 
periclustron again at about 140 Myr. Simulations A_12 to 
A_17 approach apoclustron at an age of about ~ 40 Myr, 
afther which their separations decreases again. These clus- 
ters can survive for a long time as binaries as they spiral- 
in only rather slowly. Once stellar evolution becomes less 
important (at a turn-off mass of ~ 8 M0 i.e; after about 
40 Myr) the dynamical redistribution of angular momentum 
becomes dominant again in driving the orbital evolution of 
the binary. 

In simulation model A_10 (lower solid curve in 
the stellar mass loss cannot compete with the loss of orbital 
angular momentum. The two clusters merge as soon as the 
gyro-dynamical instability sets in, which happens as soon as 
the gyration radius of the two clusters exc eeds the orbital 
separation (see Eq.l7 of lMakino et alJil99j) ). This happens 
at about 50 Myr. Stellar evolution in this case delays the 



merger by about a factor of two, but cannot prevent it. The 
gyration radius of a self gravitating system of stars can to 
first order be approximated with rvir. 

Clusters with an initial separation between about 12 
and 18 pc, have the best change to survive to an age of 
~ 100 Myr. In this range of initial separations the orbital 
evolution is strongly affected by stellar mass loss, but also by 
the redistribution of angular momentum, making this regime 
the most interesting to study numerically. If the orbital sep- 
aration exceeds about 20-30 pc the background tidal field 
starts to perturb the orbital evolution of the cluster binary. 
In addition, clusters with such wide separations are also vul- 
nerable to being ionized by passing giant molecular clouds, 
making such wide binaries unlikely to survive for an exten- 
sive period of time. 

In the range in orbital separation between about 12 pc 
and 18 pc both clusters fill their respective Roche lobes 
throughout the simulation, and stars stream through the 
first Lagrangian point from the secondary cluster to the pri- 
mary, and vice versa. Note that few stars are lost through 
the other Lagrangian points and we find no evidence for tidal 
tails in any of our simulations, contrary to what is observed 
in so me apparent binary clusters in the LMC jLeon et alJ 
I1999I) . 

The number of stars transferred from one cluster to the 
other is generally small. By the end of the simulations a to- 
tal of 238 ± 42 stars that are originally born in the primary 
cluster end up as members of the secondary cluster, whereas 
517 ± 302 stars from the secondary cluster are found in the 
primary cluster (the error here is estimated from the dis- 
persion between the different runs). The mass ratio of the 
binary remained roughly constant with time, as is illustrated 
in fig.|Hl In this figure we present the evolution of the mass 
ratio for simulation model A_20 with various initial mass 
ratios. Also in the other simulations the mass ratio hardly 
changes with time (see Tab.|5J. The simulation without stel- 
lar mass loss (dotted curve in fig.|HJ experiences the largest 
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Figure 8. Evolution of the mass ratio {q = m/M) normalized 
to the initial mass ratio for simulation model A_20. The solid and 
dotted curves are for an initial mass ratio qi = 0.167. For the solid 
curve stellar mass loss was taken into account, whereas the dotted 
curve presents the simulation in which stellar evolution was ig- 
nored. The dashed and dash-3-dotted curves are for qi = 0.33 and 
qi = 0.67, respectively. The evolution of the mass ratio for sim- 
ulation model A_20 with gi = 1 is statistically indistinguishable 
for the dash-3-dotted curve. 




t [Myr] 

Figure 9. Evolution of the separation for simulation models A at 
an initial separation of 20 pc for various choice of the initial mass 
ratio. The solid and dotted curves are for the simulation models 
A_20 with a mass ratio of 0.167 while including stellar mass loss 
(solid curve) and without stellar mass loss (dotted curve). The 
dashed and dash-3-dotted curves are for a mass ratio of 0.33 and 
0.66, respectively, while taking the total number of stars the same. 



variation in mass ratio, as in this model stars effectively 
stream from the secondary to the deeper potent ial well of 
the primary cluster (see also iMakino et al] (^^)). 

For completeness, in fig.|5]we present the evolution of 
the orbital separation for some of the simulation models pre- 
sented in fig.|Sl 

The behavior of the cluster binary is substantially af- 
fected by stellar mass loss, and the lower limit of the ini- 
tial mass function is therefore an important parameter (see 
i) l2.2ll . Mass loss in the early stage is governed by the choise 



of m_ and the slope of the IMF. In the next section we will 
further explore the range in initial cluster separations for 
which stable binaries exist by varying the choise of m_. 

3.3 The evolution of simulation models B, C and 
D 

Mass loss in the early evolution of binary star clusters is of 
major importance (see figs.|2101Eland|SJ. Most mass is lost 
by the evolving stellar population and a different initial mass 
function may therefore have dramatic consequences to the 
orbital evolution of the binary. One of the effects of mass loss 
is that the simulated clusters in model A were substantially 
less massive than observed. To compensate for this we in- 
crease the initial mass of the primary cluster from 27500 Mq 
to 40800 M0 . Since the mass ratio for the more realistic mod- 
els remains roughly constant with time (see fig.|HJ we keep 
the same mass ratio of g = 0.167. 

Much of the mass loss behavior in our simulation mod- 
els is dominated by the choise of the lower limit to the 
initial mass function. Our adopted lower limit of m_ = 
1 M (7) can be cons i dered on the high side ( see however 
SmT th fc GallaeheJ ll200ll) : iFiger et al.l J2002h : IStolte et alJ 
1I2OO5I) '). We study the effect of reducing m_ while increas- 
ing the number of stars to keep the total mass constant, but 
keeping all other parameters as much as possible the same. 
For a minimum mass m_ = 0.63 Mq we requires a total 
of 22800 stars, and for m_ = 0.50 M0 the total number of 
stars in the simulation increases to 29000. For practical rea- 
sons we do not simulate larger clusters, tho ugh in pr i nciple 
it is possible to redo these simulations with a lKroupal (I2OO j) 
initial mass fucntion with the hydrogen burning limit as a 
minimum mass. This would increase the number of stars in 
our simulation to ~ 70000, which is beyond the scope of the 
current paper. 

The evolution of the mass of the primary cluster for 
models C_20 and D_20 are presented in fig.|Sl and both 
match the observed mass of NGC 2136 and NGC 2137. We 
present the evolution of the orbital separation for these runs 
in fig. llOl For comparison we also show models A_20 with 
and without stellar mass loss. Due to the larger proportion 
of low mass stars in simulation model C_20 and D_20 stellar 
mass loss is less important than in model A_20. As a conse- 
quence the orbit (and clusters) expans less dramatically. In 
addition, the increased number of stars, and therefore the 
number density in the cluster, allows for angular momentum 
to be carried away by escapers more effectively. The result 
is a less pronounced increase in the orbital separation com- 
pared to model A_20 (see S I3.2|I . The induced orbital eccen- 
tricity is also smaller. For comparison with the observations 
we over-plot the projected distance between NGC 2136 and 
NGC 2137 in Fig.CSl 

Over-plotted with the thin solid and dash-3-dotted 
curves in fig. llOl are the result of orbital expansion by solving 
the two-body problem which includes stellar mass loss (see 

sum 



The lower (thin) solid curve gives the orbial evolu- 
tion of the 2-body results of model A_20, whereas the thick 
curve gives the results of the full A''-body simulation of the 
same model. We also show the 2-body and the full A'^-body 
results for model C_20 (see the thin and thick dash-3-dotted 
curves, respectively). To illustrate the effect of the redistri- 
bution of angular momentum we also show the result of the 
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Figure 10. Evolution of the separation for simulation models 
at an initial separation of 20 pc. The top (thick solid) curve is for 
simulation model A_20, the middle (thick dash-3-dotted) curve is 
for model C_20. The bottom (dotted) curve gives the evolution 
of the orbital separation for model A_20 but without stellar mass 
loss. The thin solid and dash-3-dotted curves give the results of 
our 2-body simulations of ^ I2.2.2l with identical initial conditions 
as for models A_20 and C_20 respectively. The mass loss for these 
models are presented in fig.|S| The observed projected separation 
between NGC2136 and NGC2137 is printed to the right as a lower 
limit. 

full A'^-body simulation of model A_20 but without stellar 
mass loss (dotted curve). 

One expects naively that the 2-body integration should 
systematically overestimate the orbital expansion, but for 
model A_20, this is clearly not the here the 2- 

body model falls short by about 40% with respect to the 
A'^-body simulation (thick solid curve). The orbital expan- 
sion of simulation model C_20, however, is quite well repro- 
duced with the 2-body solution. Remember here, however, 
that the semi-analytic models are rather simplistic, as the 
2-body problem is solved without accouting for the drag and 
friction between the clusters and, ignores the possibility of 
the redistribution of rotational angular moment of the clus- 
ters and the orbital angular momentum, whereas Eq.|3gives 
only an orbital average result. 

The best match between the simulations and the ob- 
served clusters is obtained with an initial separation of about 
20 pc with an initial mass function that extends to lower 
masses (0.5-0.6 Mq). For these models, also the total mass 
and the mass ratio of the clusters at an age of 100 Myr 
matches the observations quite satisfactorily. For model 
C_15 the orbital separation at an age of 100 Myr falls short of 
the observed separation. We therefore conclude that models 
C_20 and D_20 most satisfactorily reproduce the character- 
istics of the oberved cluster pair NGC 2136/NGC2137. 

3.4 Evolution of the cluster structure 

In fieure tTTl we present the evolution of the Lagrangian radii 
for simulation model C_20. Stellar mass loss has a dramatic 
effect on the evolution of the internal structure of the clus- 
ters. In particular the secondary cluster (thin lines in fig. llljl 
experiences a dramatic expansion in the first several 10 Myr 
due to stellar mass loss and internal dynamical evolution. 



t [Myr] 

Figure 11. Evolution of the Lagrangian radii for simulation 
models C_20. The thick curves are for the primary cluster whereas 
the thin lines are for the secondary cluster. From bottom to top 
the curves are core radius (dotted curves), 25% (dashes), 50% 
(solids) and 75% (upper dashes) Lagrangian radii. The top thick 
and thin dotted curves give the tidal radius for the primary and 
secondary cluster, respectively. Note that we smoothed the data 
for the secondary cluster over a co-moving window using 4 points 
spanning about 6 Myr. 

Core collapse in this cluster occurs at about 2 Myr, for 
the primary cluster a shallow core collapse is reached at 
about 70 Myr. Due to the induced rotation by their mutual 
tidal coupling, the core coUaps in both clusters may happen 
somewhat earlier than naively expected for isolated clus- 
ters. Rotation can initiate a gravogyr o catastrophe through 
the transport of angular mom entum jAkivama fc Sugimotd 
119891: ISpurzem fc Einsellll99g^ . 

The post collapse expansion of the secondary cluster 
exceeds the adiabatic expansion driven by stellar mass loss. 
Since the secondary cluster expands more quickly than the 
orbital separation it continues to over-fill its Roche lobe and 
transfer mass to the primary cluster, whereas the primary 
cluster detaches from its Roche lobe. As a result the orbital 
separation between 10 Myr and 40 Myr increases even more 
dramatically than naively expected from the rapid mass loss 
by stellar evolution alone (see fig. llOH . It is interesting to note 
that the evolution of a binary cluster is governed by the more 
rapid dynamical evolution of the secondary cluster, contrary 
to stellar binaries in which it is generally the primary star 
which evolves first. 



4 DISCUSSION 

We studied the evolution of bound pairs of star clusters by 
means of direct N-body simulations while including the ef- 
fects of stellar mass loss. The initial parameters selected 
for our study are loosely based on the binary star cluster 
NGC 2136 and NGC 2137 in the large Magellanic cloud. 

When starting with an initial cluster mass of about 
40800 M© and 6800 Mq in a circular orbit with a separation 
of 15-20 pc our simulations at about 100 Myr are consistent 
with the masses and orbital separation of the observed clus- 
ter pair. The observations are reproduced most satisfactorily 
when we adopt a Salpeter mass function with a lower mass 
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limit of 0.5-0.6 Mq, but lower mass cut-offs may also pro- 
duce a satisfactory comparison with the observations. 

A word of caution, however, is well placed here, as 
we did not cover the entire parameters space in excruci- 
ating detail. The initial density profile, not varied in our 
study, may have a profound effect on the evolution of the 
binary cluster. The main aim of this paper is therefore not 
to constrain the initial conditions for NGC 2136/NGC 2137, 
but rather to obtain a better understanding of the gen- 
eral evolution of binary star clusters. Earlier numerical 
studies of binary star clusters ignored the stellar mass 
spectrum, s tellar evolution and used approximate N-body 
techniques jS ugimoto fc MakinolllQSgl: iMakino et"al]|l99lt 
Ide Oliveira et al...l998i) . In our simulations all these effects 
are included in a self consistent fashion. 

We are pleasantly surprised by the rich dynamics em- 
bedded in the evolution of binary star clusters. In particular 
the finding that stellar mass loss in the early evolution of 
the cluster is rapid enough to be considered a shock for the 
orbital elements of the binary cluster came somewhat unex- 
pected. The consequence is that the orbits of star clusters 
with an age 10 Myr are not circular, but small eccentrici- 
ties e ^ 0.2 are induced upon the orbit, even if the two clus- 
ters are in Roche-lobe contact. The evolution of the distance 
between the two clusters is than a consequence of the orbital 
evolution (adiabatic/impulsive expansion) and the fact that 
cluster pair aproaches apoclustron. Note, however, that all 
our models started with circular orbits. This choise was to 
limit parameters space as introducing an initial eccentricity 
also requires the choise of an eccentric anomaly. 

Our initial clusters are not rotating. Since the redistri- 
bution of angular momentum turns out to be a mayor effect 
in the evolution of the cluster binary, future simulation may 
study the evolution of star cluster binaries by taking this 
extra parameter into account. 

Another interesting aspect of our study that the less 
massive cluster is more likely to overfill its Roche lobe. 
The secondary cluster experiences core collapse in an earlier 
stage than the primary clusters. The resulting expansion of 
the post collapsed secondary cluster subsequently drives the 
transfer of mass through the first Lagrangian point to the 
primary cluster. Mass transfer proceeds from the less mas- 
sive to the more massive component but the total number 
of stars transferred is relatively small (1-3%). The orbital 
evolution is dominated by stellar mass loss and by redistri- 
bution of angular momentum and by escaping stars 

During the exchange of stars from one cluster to the 
other hardly any stars are lost through the second and third 
Lagrangian points. Stars do escape the clusters through the 
first Lagrangian point and isotropically. The isotropic es- 
capers are mainly neutron stars which often receive high 
kick velocities in a supernova explosion. We therefore see 
no indication that tidal tails form in binary clusters. The 
origin of the tidal tails in several o bserved multiple clus- 
ters in the small Magell anic cloud bv lLeon et alj (^9^ and 
Ide Oliveira et al.l ll2000l) is therefore unlikely to be caused 
by the binarity of these clusters. 

Binary star clusters in the LMC with a separation 
exceeding about 20 pc are vulnerable to ionization by gi- 
ant molecular clouds. Our initial conditions approach this 
boundary. Such widely separated clusters can survive be- 
cause the encounter rate between clusters in the LMC is 



rather small, making it unlikely that a cluster binary is 
ionized within a few times 10 Myr. However, as we demon- 
strated with our simulations, the orbital separation may in- 
crease considerably, making such clusters much more vul- 
nerable to close (ionizing) encounters with other clusters or 
giant molecular clouds. Ionization by the background tidal 
field of the LMC is also possible. We did not take these ef- 
fects into account in our simulations, but the consequences 
may be profound. One consequence is that there are (at 
least) two reasons why binary clusters in older populations 
are rare: close cluster pairs tend to merge and wide cluster 
pairs are likely to be ionized by the background tidal field or 
by passing Molecular clouds. Of course, we cannot exclude 
that binary clusters are simply rarely formed. 

If the distance between the two clusters is initially suf- 
ficiently large that the redistribution of angular momentum 
and escapers has little effect on the orbital evolution, semi- 
analytic calculations about the orbital evolution of the bi- 
nary cluster produce satisfactory results. As soon as dy- 
namical effects start to become important (at separations 
^ 15 pc) these approximations break down. 
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